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SECTION  I 
INTRODUCTION 

The  Air  Force  Weapons  Laboratory  (AFWL)  has  been  tasked  to  calculate  the 
blast  loading  on  a shelter  due  to  a burst.  For  this  problem  a shock  front  will 
be  needed  which  is  stable  as  it  propagates  through  the  air,  yet  is  sharply  de- 
fined (spread  over  as  few  zones  as  possible).  To  provide  the  stability,  we 
have  used  the  well-known  technique  of  adding  a pseudo-viscous  pressure,  or 
artificial  viscosity,  into  the  hydrodynamic  equations  (ref.  1).  Stronger 
shocks  require  larger  viscosity  coefficients  than  weaker  shocks  to  maintain 
stability,  but  using  the  larger  coefficients  in  the  weaker  shocks  smears 
the  shock  front  over  many  zones.  This  is  unacceptable  for  the  shelter  calcu- 
lations since  the  peak  overpressure  on  the  shelter  as  a function  of  time  is 
required.  As  the  calculations  will  be  made  for  many  pressure  levels,  the  AFWL 
has  developed  a formula  for  the  linear  viscosity  coefficient  as  a function  of 
pressure  ratio. 

A preliminary  problem  set  was  devised  in  order  to  generate  and  evaluate 
the  incident  waveform  over  a wide  range  of  pressures.  In  order  to  generate 
the  incident  square  wave,  the  left  boundary  conditions  for  the  HULL  (ref.  2) 

2-0  code  were  modified  to  include  a subroutine  which  used  Rankine-Hugoniot  re- 
lations to  specify  the  density,  horizontal  particle  velocity  and  internal 
energy  of  air  for  any  incident  pressure  level.  For  the  pressure  ranges  in 
this  report,  y (the  ratio  of  specific  heats)  was  a constant  = 1.4. 

The  ambient  values  used  were: 


pressure 

Po 

1 .00156E6 

dynes/cm^  = 14.53  psi 

density 

*^0 

1.22499E-3 

gm/cm^ 

sound  speed 

3.38327E4 

cm/ sec 

energy 

^0 

2.044E9 

ergs/gm 

The  parameters  used  to  specify  the  incident  waveform  at  various  test  pressures 
were: 
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P 

■(P.si) 

2 

20 

200 


RHO 

(qtn/cm^) 
1.343447E-3 
2.2316484E-3 
5. 285441 2E-3 


I 

(erqs/qm) 
2.1208678E9 
2.6666835E9 
6. 995801 9E9 


U 

(cm/sec) 
3.1466528E3 
2. 2533251 E4 
9.2991486E4 


where  P is  the  incident  overpressure,  RHO,  I,  and  U are  the  density,  specific  internal 
energy,  and  velocity  of  the  material  behind  the  shock  front,  respectively.  The 
variable  P is  the  desired  incident  overpressure  and  is  specified  in  subroutine 
LBOUND.  Variables  RHO,  I,  and  U are  input  variables  for  a package  of  air. 

These  input  variables  are  determined  from  Rankine-Hugoniot  relations  where  P^ 
is  the  incident  overpressure,  cs  is  the  ambient  sound  speed,  and  y = 1-4: 


sc 


p = p, 


r 7 - ® Pso/Poi 
•L  7 Pso/Po 


5 P 


U = 


so 


7 P. 


cs 


0 + 6 P^„/7 


and 


I = 


P + P 
0 so 


A two-dimensional  rectangular  universe  was  generated  with  x = 200  cm  and 
y = 10  cm,  forming  a matrix  of  100  x 6 zones.  The  timeframe  of  interest  was 
10  uS  to  5 ms.  In  this  way  both  the  incident  and  reflected  waveforms  could  be 
observed,  giving  information  concerning  the  stability  and  shape  of  the  waveform 
at  different  pressures. 
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SECTION  II 

TWO  PSI  CALCULATIONS 


DEFAULT  WAVEFORM 

For  the  first  problem  (19.8010)  the  incident  pressure  was  the  default  value, 
2 psi.  This  problem  was  run  to  check  the  left  boundary  conditions  and  to  eval- 
uate the  default  waveform  (no  artifical  viscosity).  The  HULL  code  adjusts  its 
time  step  using  the  Courant  condition  which  says  simply  that  no  signal  may 
propagate  a distance  greater  than  the  smallest  zone  size  in  one  time  step.  For 
further  control,  the  user  may  set  a variable  STABF  (stability  factor)  which 
multiplies  the  Courant-determined  time  step  to  calculate  the  final  time  step 
used  by  the  code.  For  the  default  waveform,  the  stability  factor  was  0.5.  The 
left  boundary  conditions  (where  the  waveform  is  generated)  were  seen  to  be  quite 
stable.  The  incident  waveform  oscillated  somewhat  at  early  times  due  to  the 
finite  differencing,  but  when  the  wavefront  hit  the  wall  (at  2.5  msec),  much  of 
the  oscillation  was  smoothed  out  and  the  pressure  overshoot  was  about  0.7%  (fig, 
1).  The  leading  edge  was  smeared  over  about  ten  zones.  The  reflected  wave 
(4.2  psi  overpressure)  was  smeared  over  approximately  15  zones,  with  an  over- 
shoot of  about  0.9%  at  4.8  msec  (fig.  2). 

THE  EFFECT  OF  THE  STABILITY  FACTOR 

A second  problem  (19.8011)  was  run  with  the  same  input  conditions  as  the  de- 
fault waveform  except  that  the  stability  factor  was  decreased  by  an  order  of 
magnitude  to  0.05,  also  with  no  artifical  viscosity.  The  wave  propagated  at 
about  the  same  velocity  as  the  default  waveform,  but  required  over  six  times 
the  number  of  cycles  in  order  to  reach  the  wall.  The  early  oscillations  did 
not  smooth  out  as  the  waveform  moved  down  the  tube.  Instead,  the  waveform  be- 
came increasingly  unstable.  The  pressure  overshoot  for  the  incident  wave  at 
2.5  msec  was  about  2.5%  (fig.  3)  and  by  4.8  msec  the  reflected  wave  overshoot 
was  about  3.8%  (fig.  4).  The  incident  wave  (2.5  msec)  was  smeared  over  about 
12  zones.  We  conclude  that  the  default  stability  factor  (0.5)  is  more  suitable 
for  these  problems,  giving  a better  waveform  at  reduced  cost  (because  fewer 
cycles  are  required  for  propagation). 
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DEFAULT  ARTIFICIAL  VISCOSITY 

The  third  problem  (19.8012)  was  run  with  the  artificial  viscosity  option 
turned  on.  This  option  provides  a pseudoviscous  term 

q = pAu  (-  Cl  cs  + Cq  au) 


where  cs  is  the  sound  speed,  au  is  the  velocity  gradient,  and  the  constants  Ci 
and  Cg  are  the  coefficients  of  linear  and  quadratic  viscosity,  respectively. 

This  term  is  added  to  the  pressure  at  each  cycle  (ref.  1).  The  default  coef- 
ficients are  Ci  = 0.5  and  C^  = 0.  The  resultant  waveform  (problem  19.8012)  was 
found  to  propagate  slightly  faster  than  the  waveform  with  no  viscosity  (fig.  5). 

By  the  time  the  wavefront  had  reached  the  wall  (2.5  msec)  the  early  oscillations  j 

were  gone,  with  no  overshoot  or  undershoot  noticed.  The  only  difficulty  was 

that  the  wavefront  was  rather  rounded  and  spread  over  quite  a few  zones— 15 

zones  for  the  incident  wave  at  2.5  msec  and  > 30  zones  for  the  reflected  wave 

at  4.8  msec  (fig.  6).  Though  pretty,  the  waveform  was  not  really  close  to  a 

i 

square  wave.  We  decided  to  vary  the  viscosity  coefficients  in  an  attempt  to  i 

improve  the  waveform.  | 

VARIATIONS  IN  THE  VISCOSITY  COEFFICIENTS  | 

f 

Three  problems  were  then  run  with  different  viscosity  coefficients.  Prob-  I 

lem  19.8010A  was  run  with  a small  linear  coefficient  (Ci  = 0.05  and  C^  = 0.), 

The  problem  19.8011A  was  run  with  a large  quadratic  coefficient  (Cj  = 0.  and  N 

= 0.5),  and  problem  19.8012A  was  run  with  a moderate  linear  coefficient  h 

0 j 

(Cl  = 0.2  and  C^  = 0.).  Adding  the  quadratic  term  in  19.8011A  improved  the 

original  (default)  waveform  only  slightly  (figs.  7 and  8),  reducing  the  re-  jj 

i • 

fleeted  overpressure  overshoot  to  about  0.7%  while  being  smeared  about  the  same  i| 

(10  zones).  Adding  a small  linear  term  (19.8010A)  further  improved  the  orig- 

inal  waveform,  reducing  the  reflected  overshoot  to  0.4%  while  maintaining  the  j| 

same  smear  (figs.  9 and  10).  Adding  a moderate  linear  term  (19.8012A)  im-  j- 

proved  the  original  waveform  by  eliminating  the  overshoot,  but  at  the  expense  j 

of  increasing  the  smear  to  twelve  zones  (figs.  11  and  12).  We  consider  the 
small  linear  term  to  be  the  better  alternative. 
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AN  ARTIFICIAL  VISCOSITY  FUNCTION 

An  artifical  viscosity  function  was  developed  to  provide  an  appropriate 
artificial  viscosity  coefficient  for  any  incident  overpressure  between  2 psi 
and  about  1200  psi.  This  formulation  is  an  empirical  fit  applied  to  SAP  (ref. 
3,  a 1-D  Lagrangian  code)  to  develop  the  1-KT  Standard  (ref.  4).  A slowly- 
varying  function  of  pressure  was  needed  which  approaches  an  upper  limit  as 
pressure  increases.  Density  as  a function  of  pressure  satisfies  these  cri- 
teria. The  linear  coefficient  of  viscosity  was  formulated  to  be 

Cl  = /(O.l)  Rhm  + 0.015 


where 


Ohm  - 1 

Rhm  - — - 1 . 

o 


> p. 


and  p„^  is  the  maximum  density  found  in  the  grid  while  p is  the  ambient  dens- 
rndx  o 

ity.  The  quadratic  coefficient  is  given  by 

'o  = ° 

With  this  artificial  viscosity  function  (AVF)  implemented  (problem  19.8012B), 
Cl  was  approximately  0.113  (incident)  and  0.156  (reflected)  (figs.  13  and  14). 
With  these  coefficients,  the  2 psi  waveform  was  comparable  to  the  best  of  the 
preceding  runs,  problem  19.8010A.  On  the  incident  wave  the  overshoot  was  only 
0.4%  and  the  shock  was  smeared  over  10  zones  at  2.5  ms.  At  4.8  ms  the  over- 
shoot was  less  than  0.2%,  while  the  shock  smeared  15  zones. 
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SECTION  III 

TWENTY  PS  I CALCULATION 

WITH  THE  ARTIFICIAL  VISCOSITY  FUNCTION 

Only  a single  problem  (19.8013)  was  run  with  a 20  psi  incident  waveform,  to 
check  the  AVF  at  20  psi  and  at  about  60  psi  (the  reflected  wave).  For  the  inci 
dent  wave,  Ci  = 0.302  (incident)  and  0.474  (reflected).  The  incident  waveform 
at  1.8  msec  (just  before  impact)  shows  no  overshoot  and  is  spread  over  only 
nine  zones  (fig.  15).  The  reflected  waveform  at  4.2  msec  shows  a very  slight 
overshoot  and  an  acceptable  spread  of  10  zones  (fig.  16). 
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Figure  16.  AFWL  - 20  psi  Square  Wave  on  Wall 
Time  - 4.2  msec,  Problem  19.8013. 
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SECTION  IV 

TWO  HUNDRED  PSI  CALCULATIONS  i 

! 

DEFAULT  WAVEFORM  j 

The  first  200  psi  problem  (19.8014)  was  run  to  evaluate  the  default  wave- 
form (stability  factor  of  0.5,  with  no  artificial  viscosity).  At  200  psi  the 
default  waveform  is  nearly  perfect,  having  virtually  no  overshoot  and  a spread 
of  four  or  five  zones  just  before  hitting  the  wall  at  700  ysec  (fig.  17).  The  ! 

reflected  wave  at  2.5  msec  (1200  psi),  however,  was  quite  poorly  formed,  with  I 

an  overshoot  of  about  25%  (fig.  18).  ■ 

THE  EFFECT  OF  DEFAULT  VISCOSITY  | 

P 

The  second  200  psi  problem  (19.8014A)  was  run  with  the  default  viscosity  : 

coefficients  (Cj  = 0.5  and  C^  = 0).  The  incident  waveform  at  700  usec  had  no  j- 

overshoot  but  was  spread  over  about  eight  zones,  which  was  acceptable  (fig.  19).  ! 

The  reflected  wave  at  2.5  msec  was  considerably  improved  over  the  default  200  |'; 

psi  problem,  having  an  overshoot  of  about  1.8%  and  a spread  of  six  zones  (fig.  I 

20).  II 

i i 

WITH  THE  ARTIFICIAL  VISCOSITY  FUNCTION  1 I 

I i 

Finally,  a third  200  psi  problem  (19,80148)  was  run  with  the  AVF  implemented.  ji 

For  this  problem,  Cj  = 0.591  (incident)  and  1.106  (reflected).  The  incident  j| 

(200  psi)  waveform  at  700  usec  showed  no  overshoot,  a slightly  rounded  top  lead-  j 

ing  edge,  and  an  acceptable  spread  of  about  eight  zones  (fig.  21).  The  re-  1 

fleeted  waveform  (at  about  1200  psi)  at  2.5  msec  showed  virtually  no  overshoot,  ^ 

a slightly  rounded  bottom  leading  edge,  and  a spread  of  about  10  zones  (fig.  22).  ] 


27 


AFWL-TR-77-53 


AFWL-TR-77-53 


Pf>CS«iK  r«CRi;0*4Tn»  nlSTCOPW^ 

CPJSS  VCTiCf4  ;j?  PilHur-t  ‘».rC-05  KM  (j: 


I 20  hO  60  ao  ICO 


Mtn  I i.oeM«» 

mt : 

Oai  : ?.axi>CQ 


H(3m?nNirv  vClCCIM  hcth/jhiol 
HilMl/ONliic  C»‘0.»S  iiCHOfl  K JlP  Jt  •*.  ri  O'j  i j: 

fi  ^?5”  ~m5  60  liJ  icb 


S *•“  1 


:j 

•*•1  ; I.  r. *•;« 

ni  ; ; Ju4>x 


•*  •!  •!  i.>  \.%  .M  Kt 


0.0  .1  • 


i.>  l.«  .'.I 


Ul^fSITT  i-iCHI/UNiat.  mIMOOJ'M* 

Mmi/ONtiX  rnoss  'tenon  ST  (XniuOt  l.lTt-OS  fM  ij: 

1 Hi  ^0  60  eo  ioo 

***  t l.;>4-9] 

' Oil  ; 


CNt^^'iT  ‘Cni/Cirfle 

MOHUflurrtu  cws  s£cn;:M  st  s.i:rL:{  >4.ri-us  "m  •. 

fi  n ho  S c5  Tco 


i ;.Sv«4-'i 
Ut  : } 


*•  ■>  ••  .•  >.i  *.v  i.«  /.I 


.1  .»  .1  i.V  ».%  i.«  .Ji  ».♦ 


Figure  21.  AFWL  - 200  psi  Square  Wave  on  Wall  - Viscosity  Function 
Time  - 700  usec,  Problem  19.30148. 
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SECTION  V 


REFLECTION  FACTORS 

One  method  of  validating  the  results  of  these  calculations  is  to  check  the 
value  of  the  peak  reflected  overpressure  for  each  problem  against  the  peak  for 
an  ideal  gas  (y  = 1.4).  Using  Rankine-Hugoniot  relations,  the  formula  for  peak 
reflected  overpressure  is 


aP  = 2aP  { ^'*‘^4P/P  \ 

\ 7+aP/P  ) 

where  aP  = P,.-P,  and  aP  = P -P,  (P.  = incident  pressure,  P = reflected  pres- 
ia  rrai  r 

sure,  and  P,  = ambient  pressure). 

d 

The  following  table  summarizes  the  results  for  those  calculations  which 


used  the  AVF: 

AP 

aP  (actual 

r 

aP^( theory) 

Problem 

(dynes/cm^) 

(dynes/cm2) 

(dynes/cm^) 

19.8012B 

1.381E5 

2.924E5 

2.923E5 

19.8013 

1.379E6 

4.132E6 

4.131E6 

19.8014B 

1.379E7 

8.234E7 

8.249E7 

There  is  a 0.2%  error  in  the  200  psi  problem  and  better  agreement  with  the 
others.  This  clearly  suggests  that  the  HULL  code  is  modeling  the  head-on 
shock  properly. 


SECTION  VI 


CONCLUSIONS 

The  HULL  code  can  provide  a stable,  sharply-defined  shock  for  any  over- 
pressure level  from  2 psi  to  1200  psi  using  the  artificial  viscosity  function 
presented  in  this  report.  Both  the  incident  and  the  reflected  waveforms  are 
stable  yet  defined  sharply  enough  so  that  valid  pressure  levels  for  blast  en- 
vironments may  be  calculated. 
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